A three-state cyclic voter model extended with Potts energy 
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The cyclically dominated voter model on a square is extended by taking into consideration the variation 
of Potts energy during the nearest neighbor invasions. We have investigated the effect of surface tension 
on the self-organizing patterns maintained by the cyclic invasions. A geometrical analysis is also developed 
to study the three-color patterns. These investigations indicate clearly that in the "voter model" limit the 
pattern evolution is governed by the loop creation due to the overhanging during the interfacial roughening. 
Conversely, in the presence of surface tension the evolution is governed by spiral formation whose geometrical 
parameters depend on the strength of cyclic dominance. 

PACS numbers: 02.50.-r, 05.50.+q, 87.23. Cc 



I. INTRODUCTION 

In systems with several species (particles, opinions, 
etc.) the cyclic invasion processes can maintain a self- 
organizing domain structure. Although this phenomenon 
is investigated extensively in different areas, such as 
the chemical reactions on crystal surfaces biolog- 
ical (Lotka-Volterra) systems Rock-Scissors-Paper 
(RSP) games in evolutionary game theories 0, cyclically 
dominated voter models |g|-[lC|], the mechanism sustain- 
ing the polydomain patterns is not well understood. At 
the same time, this type of spatial self-organizations is 
believed to play crucial role in the biological evolution 
|llf and it can provide protection for the participiants 
against some external invadors Jl2|;[l3| • 

One of the simplest model exhibiting variations in 
a self-organizing domain structure was introduced by 
Tainaka and Itoh J§|,[)| . In this cyclically dominated voter 
model three states {A, B, and C) are permitted on the 
sites of a square lattice. The system evolution is governed 
by the iteration of invasions between two randomly cho- 
sen nearest neighbors. The cyclically symmetric invasion 
rates are characterized by a single control parameter P. 
For example, the pair AB transforms into AA with a 
probabilty 1/2 + P or to BB with probability 1/2 - P 
(0 < P < 1/2). For P — this system is equivalent 
to the voter model p5|JTf^ | exhibiting growing domains 
whose correlation length is proportional to yt if the sys- 
tem is started from a random initial state. Contrary, for 
P = 1/2 A beats B beats C beats A as it happens in 
the RSP game. In this former case the above rules sus- 
tain a self-organizing, three-color polydomain structure 
in which the linear domain size can be characterized by a 
correlation length of £ ~ 2.5 (measured in lattice unite). 
Tainaka and Itoh have shown that the typical domain 
size diverges when P — > Their numerical analysis 
is focused on the density of vortices which are defined by 
those points of a three-color map where the three states 
(domain walls) meet. In fact, one can distinguish vortices 
and antivortices (rotating clockwise and anti-clockwise) 



and they are created and annihilated in pairs during the 
evolution of domain structure. According to the early 
Monte Carlo (MC) simulations the vortex density can 
be approximated by a power law behavior in the limit 
P->0@. 

In the above system the domain walls separating two 
homogeneous domains are very irregular. This irregu- 
larity prevents the observation of expected spirals to be 
formed by the rotating vortex arms for smooth interfaces 
|p|,p^| . In the absence of cyclic invasion (and related vor- 
tex rotation) some features of the rough interfaces are 
already studied by several authors considering the two- 
state voter model ^5|-^o| . In the presence of cyclic inva- 
sion, however, the topological and geometrical features 
of domain structure are not yet investigated rigorously 
even for smooth interfaces. Very recently some geometri- 
cal features and drift of a single spiral wave are studied by 
using a continuous reaction-diffusion model |2l],|22| . Now 
we report a model which exhibits a transition between 
different self-organizing patterns involving those where 
smooth, rotating spiral arms can be observed. 

The present work is devoted to study the effect of 
surface tension on the self-organizing domain structures 
maintained by cyclic invasions on a square lattice. For 
this purpose we have introduced a model where the near- 
est neighbor invasion is affected by an interfacial (Potts) 
energy |2^,Q whose strength is controlled by a param- 
eter K. Using MC simulations we have investigated the 
vortex density in the stationary state for different val- 
ues of P and K. In order to have a deeper and more 
quantitative insight into the domain structure we have 
developed a method to study some geometrical features 
of the interfaces. This analysis confirms the necessity of 
these types of sophisticated descriptions. 



II. THE MODEL 

We consider a square lattice where at each site x = 
(i and j are integers) there is a state variable with 
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three possible states, namely s x — 0, 1, 2. The Potts 
energy for a configuration s — {s x } is defined as 



H = — 



l S ( S X,Sy) - 1] 



(1) 



<x,y> 



where the summation runs over the nearest neighbor sites 
and 6(s,s') indicates the Kronecker's delta f^jjj^ij] . No- 
tice that the coupling constant is chosen to be energy 
unit. In the present form the Potts energy measures the 
length of equivalent interfaces (in lattice unit a = 1) sep- 
arating the three types of domains and its inverse es- 
timates the average domain radius p5[ . Evidently, in 
the threefold degenerated (homogeneous) ground state 
H = 0. 

The configuration evolves in time according to elemen- 
tary invasions between two nearest neighbor sites (x and 
y) chosen randomly. More precisely, a pair of neighbor- 
ing state variables (s x , s y ) (assuming s x ^ s y ) transforms 
into (s y ,Sy)) with a probability 



no 



where 



x, 



(Sy,S y )] = 



1 



1 + cxp (KAH + PD) 



AH = Hf- Hi 



(2) 



(3) 



is the energy difference between the final and initial 
states, and K, as an inverse temperature, controls the 
effect of Potts energy on this single site flip. The second 
term in the argument of exponential function describes 
the cyclic dominance with a strength P, where 



D = 



+ 1 if s x = (s y 
-1 if s x = (s y 



1) mod 3 , 

2) mod 3 . 



(4) 



In the case K = the present model is equivalent to those 
introduced by Tainaka and Itoh [|| . Evidently, the three- 
state voter model is reproduced if P = and K = 0. For 
K > 0, however, the interfaces become more smooth be- 
cause the present dynamics suppresses those elementary 
processes which increase the interfacial (Potts) energy. 

We have to emphasize that for P = this system ex- 
hibits domain growth independently of the value of K. 
Apparently, for P = Eq. (^) satisfies the condition of 
detailed balance at a temperature 1/K as defined for 
the kinetic Potts model (on the analogy of kinetic Ising 
model) that undergoes a symmetry breaking phase tran- 
sition when increasing K. There exists, however, a rele- 
vant difference due to the constraint of invasion dynam- 
ics. Namely, in the present model the new state at a 
given site should be equivalent to one of the neigboring 
one. This means that the changes are localized along the 
boundaries separating the homogeneous domains. 

Notice that the above rules conserve the cyclic symme- 
try among the three states. As a result in a sufficiently 
large system the three states are present with the same 
probability (1/3). For small system, however, one of the 



species can extinct due to the effect of fluctuations and 
finally the system evolves toward one of the three (homo- 
geneous) absorbing states. Henceforth our analyses will 
be restricted to the large system limit that is provided 
by choosing the system size to be much larger than any 
length characteristic to the corresponding pattern. 

The above model is investigated by MC simulations 
under periodic boundary conditions on a square lattice 
consisting of N — L x L sites. The system is started 
from a random initial state where the three states are 
present with equal probabilities. During the simulations 
we have recorded the number of vortices and antivortices 
defined above. For this purpose we have counted those 
2x2 block co nfig urations containing all the three pos- 
sible states PHlOll . After a suitable transition time we 
have determined the ave rag e vortex density as well as 
its fluctuation defined in []l Cf| . The system size is varied 
from L — 400 to 2000 to have sufficiently large number of 
vortices in the stationary states. The numerical results 
of vortex densities are summarized in a log-log plot as 
demonstrated in Fig. [[} 
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FIG. 1. Vortex densities as a function of P for K = 1 
(closed triangles) , 1/4 (open triangles) , 1/16 (closed squares), 
1/64 (open squares), and (pluses). The solid line shows the 
predicted power law behavior if K = 0. The dashed line (with 
a slope of 2.05) indicates the best power law fit for K = 1. 

If P >> max(A", 1) then the dynamics is governed 
by the deterministic RSP rule that maintain a self- 
organizing state with small domain sizes (£ ~ 2.5) as 
mentioned above. Consequently the vortex density be- 
comes independent of K for sufficiently large values of P 
as demonstrated in Fig. |l|. 

In the case K = the P-dependence of vortex density 
can be well described by a power law, namely, p v ~ P^ 
within the range 0.003 < P < 0.3. The best fit is found 
for P — 0.29(1) confirming our previous result |I(J. At 
the same time Fig. [j] demonstrates clearly that the vortex 
density is dramatically reduced when the interfacial en- 
ergy is switched on. For K = 1 the MC data can be well 
approximated by another power law with an exponent 
(3 = 2.05(9). It should be emphasized that within the 
statistical error our data is consistent with a quadratic 
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behavior. Similar behavior can be conjectured from the 
trends represented by MC data for lower K in Fig. |l|. Un- 
fortunately, we couldn't confirm this expectation by de- 
termining the leading term in the P-dependence of vortex 
density for lower K values because this analysis requires 
extremely long run time and large sytems. Just to in- 
dicate the difficulties, the determination of a data point 
at low vortex densities has required more than four-week 
run time on a fast PC. We think that further numerical 
analyses are necessary to justify (or modify) the above 
conjecture. 

Due to the long run times we could derive the vor- 
tex density fluctuations (y) with an adequate accuracy. 
As demonstrated in Fig. || the numerical data indicate 
the divergency of the vortex density fluctuation in the 
absence of interfacial energy (K = 0). Within the in- 
vestigated region, this fluctuation can be approximated 
as x ~ P -7 with 7 = 0.3(1) in good agreement with a 
previous result Q . 
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FIG. 2. Fluctuation of vortex density versus P for different 
K values denoted by the same symbols as in Fig. ^ The solid 
line indicates the fitted power law divergency in the absence 
of interfacial energy. 

In the presence of interfacial energy (K > 0) the vortex 
density fluctuation vanishes with P as indicated in Fig. |^. 
The vortex density fluctuation seems to be proportional 
to the vortex density for sufficiently low p v . Similar fea- 
tures have been found for some systems of particles and 
antiparticles performing branching annihilating random 
walks |l4|,[l(J. This is the reason why we have reinvesti- 
gated a parallel drawn between the vortex dynamics in 
present model and a system of particles and antiparticles 
as suggested in [jio). According to a simple idea the ro- 
tating vortices form spirals whose long and narrow arms 
enhance the probability of the creation of a new vortex- 
antivortex pair. The movement of vortices can be well 
approximated by a random walk on a lattice. Further- 
more, a vortex and antivortex annihilate each other when 
meeting at the same site during their random walks. The 
balance between the annihilation process and pair cre- 
ation yields an average concentration in the stationary 
state. Within this framework a simply mean-field anal- 



ysis (details are given in |T(J) predicts that a quadratic 
behavior (p v ~ P 2 ) can be reproduced if the pair creation 

is proportional to Pp\ or P 2 p v . From the view-point of 
vortex dynamics both possibilities demand a better un- 
derstanding about the relationship between the creation 
of vortex- antivortex pairs and the geometry of interfaces. 

As mentioned above the whole length of interfaces is 
equivalent to the Potts energy defined by Eq. (|B|) . During 
the simulations we have determined the expected value 
of Potts energy per sites, 

E = ^{H) (5) 

where (• • •) indicates the average over the sampling time. 
The results of our simulations are summarized in a log- 
log plot (see Fig. ||). 
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FIG. 3. Average Potts energy per sites as a function of 
P for the same values of K plotted in Fig. jl]. The solid 
and dashed lines indicate the slopes of 0.08 (K — 0) and 
1.3 (K = 1). 

At the first glance the P- and if-dependences of the 
vortex density and Potts energy seem to be very simi- 
lar. However, the detailed numerical analysis gives dif- 
ferent values for the exponents when fitting a power law 
(E = aP a ) for small P values. Namely, we have obtained 
a = 0.08(2) and 1.3(1) for K = and 1 respectively. 
The reader can easily check that the traditional scaling 
argument (predicting /3 = 2a) is not valid in the present 
cases. 

The self-organizing domain structure shows striking 
differences depending on wether the interfacial energy is 
switched on or not. For the sake of illustration two typi- 
cal patterns are shown in Figs. ^ and [| In the absence of 
interfacial energy (K — 0) the nearest neighbor invasions 
yield irregular boundaries whose overhanging results in 
small islands (loops). Their random motion, extension, 
shrinking, splitting, and fision seem to play crucial role in 
the pattern evolution as well as for the three-state voter 
model. In this former case (P = 0), however, these el- 
ementary events are not able to prevent the growth of 
domains whose characteristic linear size increases with 
time as \fi 
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FIG. 4. Typical part (50 x 50) of snapshoot in a larger 
system for K = and P = 0.01. 

Choosing a particular initial state it is already demon- 
strated that the cyclic dominance drives the vortex ro- 
tation (for P > 0) which is accompaned with spiral for- 
mation PJTc|]. In Figure ^ the rotating spirals are not 
recognizable due to the irregular interfaces. However, 
the spiral formation becomes visible when the interfacial 
roughness is reduced by the surface tension as demon- 
strated in Fig. [|. 




FIG. 5. Spiral formation sustained by the rotating vortices 
and antivortices is recognizable on a 200 x 200 part of a larger 
system for K = 1/2 and P = 0.3. 

In Figure ^| one can easily identify the vortices and 
antivortices rotating clockwise and anticlockwise respec- 
tively. This rotation creates spirals because the average 
invasion velocity is constant. We have to emphasize that 
this pattern can not be characterized by a single length 
unit (e.g. correlation length) because the main features of 
spirals (armlength, average curvature, average distance, 



etc.) depend on the model parameters. This is the reason 
why we have developed a method to study some geomet- 
rical fetures of three-color maps on a square lattice. 

III. GEOMETRICAL ANALYSES 

On a three-color, continuous, planar map the domains 
are separated by three types of smooth boundaries. Dedi- 
cated points are the vertices where three or more bound- 
aries meet. If such a map evolves smoothly then the 
appearance of vertices with more than three edges be- 
comes negligible. Thus our analysis can be restricted to 
those maps which contain only three-edge vortices and 
antivortices. However, as we show later, the qualitative 
feature of the system remains uneffected if four-leg ver- 
tices are not ignored. One can easily check that these 
vortices and antivortices are positoned alternately along 
domain boundaries [jiof . Our geometrical analysis will be 
focused on determining the average value of arclength, ro- 
tation of tangent vector, and curvature for those bound- 
aries connecting a vortex and an antivortex. Henceforth 
the rotating vertex is called vortex. 

On a square lattice the boundaries are polygons con- 
sisting of unit length parts whose tangential rotation may 
be Acj> = ±7r/2 and 0. For a given vortex edge the tangen- 
tial rotation is determined by summarizing these quan- 
tities step by step along the edge from a vortex to the 
connected antivortex. At the same time arclength is also 
obtained as the number of steps. The elementary step 
is based on the identification of the 2x2 block config- 
urations. This algorithm assumes that first we have de- 
termined the vortex positions. To reduce the statistical 
error this procedure was repeated many times during the 
simulations. 

The above algorithm is well defined if the three-color 
pattern is free of four-edge vertices. Unfortunately, the 
investigated self-organizing patterns contain undesired 
four-edge vertices (see Figs. ^ and ||). Some of them 
(involving all the three states) can be considered as a 
vortex-antivortex pair just before their annihilation or af- 
ter they creation [|l0| . The others involve only two states 
and make the paths (from vortex to antivortex) indefi- 
nite. Both types can be removed by executing an invasion 
through one of the randomly choosen four edges. Before 
the geometrical analysis all the investigated distributions 
are slightly adjusted by repeating the random invasions 
at the four-edge vertices until they vanish. Evidently, 
the effect of these modifications on the energy or vortex 
density is negligible if the typical domains are sufficiently 
large. The most relevant effect appears at K = when 
the density of vortex-antivortex pairs is approximately 
p v /6 in the whole region of P where we studied the sys- 
tem. Consequently, the P-dependence of vortex density 
remains power law after the pattern adjustment. 

After having removed the four-edge vertices the pat- 
tern becomes topologically equivalent to the continuous, 
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three-color map mentioned above. In this case we can 
distinguish two types of boundaries, namely, loops (sur- 
rounding an isolated domain) and vortex edges (starting 
at a vortex and ending at one of the connected antivor- 
tices). Using the mentioned algorithm we have deter- 
mined the average length Z av and tangential rotation </> av 
of vortex edges. 
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FIG. 6. Average arclength of vortex edges versus P for 
different K values. Symbols as in Fig. |l|. 

In Figure ^ the log-log plot of the average length of 
vortex edges shows that l av increases slowly when P is 
decreased for K = 0. Significantly faster increase can be 
observed in the presence of surface tension. The arrange- 
ment of MC data for K = 1 has inspired us to fit a power 
law, l av = aP~ x as we had done for the vortex density 
and the Potts energy per sites. Within the same region 
of P the best fit is found for A = 1.75(5). 

The above behavior is not surprising because the av- 
erage vortex distance exhibits qualitatively similar P- 
dependence. For the quantitative analysis an average 
vortex distance can be deduced from the density of vor- 
tices as d av = l/yffc. The striking difference caused 
by the introduction surface tension and cyclic invasion 
becomes visible when the pairs of data d av and l av are 
plotted on a log-log plot (see Fig. In Figure ^ the 
straight line (l av = 1.05d av ) demonstrates those set of 
domain structures that can be well characterized by a 
single length scale. For example, such a situation can 
be observed when considering the domain growth in the 
three-state Potts model below the critical temperature. 
In the absence of interfacial energy the MC data indi- 
cates significantly different relation which may be ap- 
proximated as l av ~ <i av 8 within the given region. The 
slower increase of the average length of vortex edges can 
be explained by the increasing number of those vortex- 
antivortex pairs which have two common (short) edges. 
Shuch a pair is frequently created when the moving is- 
lands meet the third type of domain during their random 
movements. 
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FIG. 7. Relation between the average vortex distance (d av ) 
and arclength of vortex edges (Z av ) for those parameters (and 
symbols) plotted in Fig. |l| Both lengths are measured in unit 
of lattice constant. 

An opposite tendency can be recognized for those cases 
where the spiral formation becomes relevant because the 
arclength of spiral arms always exceed the distance be- 
tween the corresponding vortex and antivortex. At the 
same time our data reflects that the average tangential 
rotation of a vortex edge increases with d av . This state- 
ment is supported by those data in Fig. || we obtained for 
K > 0. Particularly, for K = 1 one can observe that both 
av (see Fig. ||) and Z av (Fig. ^) increases monotonously 
when P goes to zero. Unfortunately, we were not able to 
study what happens when the average tangential rotation 
becomes larger than 2n. 
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FIG. 8. Average tangential rotation of vortex edges as a 
function of cyclic dominance (P) for different K values (sym- 
bols as in Fig. [l]). 

In Figure || the angle of tangential rotation is measured 
in unit 7r/2 which is a natural choice on a square lattice. 
One can observe that 4> av becomes practically zero in the 
P — > limit in the absence of interfacial energy. In the 
light of this result one can think that the spiral formation 
does not play a dominant role in the pattern formation 
for K = 0. At the same time, we should keep in mind 
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that the cyclic invasion (P > 0) is required to sustain the 
self-organizing domain structures, otherwise the domains 
would grow unlimited. Unfortunately, we cannot explain 
quantitatively the microscopic mechanism yielding this 
behavior. Now we can only give some additional argu- 
ments supporting the crucial role of islands as mentioned 
above. 

The total interfacial (Potts) energy can be separated 
into two parts. The first contribution comes from the 
island boundaries and the second part from the vortex 
edges. Thus the energy per sites can be written in the 
following form: 



E — Ei + 3p v l a 



(6) 



where Ei denotes the contributions of islands to the to- 
tal Potts energy E defined by Eq. (g). The second term 
indicates that the contribution of vortex edges can be 
expressed as a product of the density of (three-edge) vor- 
tices (p v ) and the average length of vortex edges (l av ). 
Using this expression we can determine the values of Ei 
from those data plotted in Figs. |], ||, and ||. The results 
of this calculation are illustrated in Fig. ||. 
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FIG. 9. The interfacial energies of islands vs. P for differ- 
ent K values. Symbols as in Fig. [j] 

For K — 1 the interfacial energy contribution of islands 
vanishes at sufficiently low values of P. This tendency 
can be observed in Fig. ^| for other K > values. The 
lower the value of K, the lower the value of P is where 
Ei becomes negligible. This means that in the stationary 
state the number of islands are reduced by the interfa- 
cial energy. This tendency can be visually checked if the 
reader compares the two snapshots shown in Figs. ^ and 

Notice, furthermore, that for K = 1 the energy contri- 
bution of islands is negligible for weak cyclic dominance 
where the P-dependence of E, p v , and l av can be well ap- 
proximated by power laws as mentioned above. The sub- 
stitution of the corresponding expressions into @ yields 
a relation between the exponents, namely a = /?— A. Our 
numerical data support this result. 

The investigation of Ei for K — shows surprising 
results. In the case of strong cyclic dominance the dom- 



inant part of interfacial energy comes from the vortex 
edges. According to our simulations the contribution of 
Ei to the total interfacial energy increases meanwhile E 
decrease (see Fig. ||when decreasing P in the investigated 
region. Since E > Ei therefore these opposite tendencies 
imply the possibility of a break point for P < 0.002. 



IV. CONCLUSIONS 

We have studied numericaly the effect of surface ten- 
sion on the self-organizing patterns maintained by cyclic 
invasions among three species on a square latice. For this 
purpose the cyclic voter model introduced by Tainaka 
and Itoh Q is extended in a way conserving the cyclic 
symmetries. In the original model the invasion between 
two (randomly chosen) nearest naighbors is not affected 
by the neighborhood. In the extended model the nearest 
neighbor invasion rate is influenced by the neighborhood 
via taking the variation of Potts energy into account. Our 
analyses are restricted to those situations (K > 0) where 
this modification favorizes those invasions which reduce 
the length of interfaces separating the domains. 

Our simulations have justified that the introduction 
of interfacial energy causes relevant changes in the ob- 
served patterns. To have a more quantitative and so- 
phisticated picture we have determined the average value 
of some geometrical features of the interfaces (e.g. ar- 
clength and tangential rotation of vortex edges). This 
method is based on the analogy to the continuous limit 
of a three-color map. By this way we could study the 
contributions of vortex edges and islands separately. In 
the light of this analysis we can distinguish three types 
of typical domain structures. 

In the deterministic limit (P >> ma.x(K, 1)) the pat- 
tern consists of small compact domains and contains 
many vortices and antivortices. In the absence of in- 
terfacial energy (K = 0) the typical domain size as well 
as the contribution of island interfacial energy increases 
when the the cyclic dominance (P) is decreased. Here the 
island creation via the interfacial roughening seems to be 
a relevant phenomenon. Conversely, in the presence of 
interfacial energy the islands vanish when P is decreased 
and vortex (spiral) rotations dominate the pattern evolu- 
tion. The transitions among these typical behaviors are 
smooth. 

Our numerical results are obtained in a limited region 
of the parameter P due to the technical difficulties ap- 
pearing for large typical domain sizes. For some cases 
(K = and I) our data can be approximated by power 
laws in a region of P. We are, however, not convinced 
that these (expected universal) behaviors remain valid 
for lower P values. For example, we don't know what 
happens when the average tangential rotation of vortex 
edges becomes significantly larger than 2tt. Deviations 
can also appear for K = at lower P values where Ei is 
expected to decrease monotonously with P. 
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The suggested geometrical analyses confirm that the 
self-organizing patterns cannot be characterized by a sin- 
gle length unit as it happens for many other systems. 
In these cases two patterns can not be transformed into 
each other by choosing a suitable length scale. In the 
presence of interfacial energy this feature is strongly re- 
lated to the appearance of spiral vortex edges whose av- 
erage tangential rotation remains unchanged during such 
a geometrical magnification. We think that this type of 
geometrical analysis joints different approaches and mod- 
els, furthermore, it motivates a theoretical effort to find 
general relations among these quantities. 
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